function [next_pos] = ForageSGSA(G, mass, current_pos)
[N,D] = size(current_pos);

epsilon = 1e-3;
    for n = 1 : N
        acc(n,:) = zeros(1,D);
        for k = 1 : numKbest
            force = zeros(1,D);
            indexKbest = sort_fitness(k, 2);
            force(1,:) = pos(indexKbest,:) - pos(n,:);
            R = norm(force(1,:));
            force(1,:) = force(1,:) * rand() * mass(indexKbest) / (R + epsilon);
            acc(n,:) = acc(n,:) + force(1,:);
        end
    end
    acc = G * acc;
    vel = vel.*rand(ps,D) + acc;
    vel = ((vel>=vlb)&(vel<=vub)).*vel+(vel<vlb).*vlb+(vel>vub).*vub;
    pos = pos + vel;
    pos = ((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*plb+(pos>pub).*pub;

end